Direct Numerical Observation of Anomalous Diffusion in Random Media 
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We present computer simulations of anomalous diffusion, (r 2 (t)) ~ at 1 ~ s , in two dimensions. The 
Monte Carlo calculations are in excellent agreement with previous renormalization group calcula- 
tions. Interestingly, use of a high quality pseudo-random number generator is necessary to observe 
the anomalous diffusion. A linear-feedback, shift-register method leads to incorrect, super-diffusive 
motion of the random walkers. 
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1. INTRODUCTION 



Certain types of physically realizable disorders cause 
anomalous sub-diffusion in two dimensions |QJ|] . In these 
unusual cases the mean square displacement of random 
walkers does not increase linearly with time, but rather 
increases sub-linearly with time, (r 2 (t)) ~ at 1 ^ 6 . The 
exponent for this scaling is continuously variable in the 
strength of disorder. In fact, the exponent can be found 
exactly for the type of singular Gaussian disorder that 
leads to sub-diffusion |3| ^] . This special disorder cor- 
responds physically to quenched, charged defects. The 
defects can either be true, electrostatic charges or topo- 
logical "charges." In either case, the quenched charges 
interact with the moving particles of interest via a long- 
range, logarithmic potential. 

Diffusion in a variety of random media has been con- 
sidered by numerical simulation. Typical disorders that 
have been investigated include non-singular potential 
disorder j|, random fluid flows fl(i|-|l2]|, fractal media 
|l3| |U , optical molasses |lrj| , and topologically disor- 
dered structures jL7). The random walk Monte Carlo 
simulation method is used in most numerical studies of 
diffusion. To date, however, there are no numerical stud- 
ies that verify the renormalization group predictions of 
anomalous diffusion in singular, two-dimensional, ran- 
dom potential fields. Preliminary results were cited in 
the review by Bouchaud and Georges |Q, but a formal 
publication did not ensue. 

This same type of quenched, "ionic" disorder subjects 
chemical reactions to transport limitations and causes 
anomalous kinetic behavior. Two dimensions is particu- 
larly interesting, since this is also the upper critical di- 
mension for bimolecular reactions |l9|-p2| , in addition to 
being the upper critical dimension for the charged disor- 
der. The kinetics of chemical reactions, of course, have 
always been a subject of interest for scientists and engi- 
neers. In most cases, the reactant diffusion is normal, and 
the effect of diffusion limitations on the rate of chemical 
reaction is easy to calculate. Anomalous diffusion, how- 
ever, leads to anomalous kinetics. In this case, the effect 
on the rate of chemical reaction is not so widely known. 
The kinetics of the reactions A + >0,.A + .B— >0, and 



A + + B ^= AB in singular disorder have been derived 

r 

analytically |2l] - Q . In these references, a field-theoretic 
treatment of anomalous kinetics was worked out, and 
renormalization group predictions were derived. These 
studies show that the reactions become transport lim- 
ited in the long time regime. At long times, where the 
diffusion is anomalous, the kinetics also becomes anoma- 
lous. 

Simulations to test these theoretical predictions would 
be of great interest. A necessary prerequisite for proceed- 
ing with these numerical studies of anomalous kinetics is 
first the ability to simulate anomalous diffusion. This 
ability requires both constructing the disorder and per- 
forming the diffusive motion in the quenched disorder. 

In this article, we present numerical observations of 
anomalous diffusion in two dimensions using the Monte 
Carlo method. In section 2 we discuss the appropriate 
form of the quenched disorder. In section 3 we intro- 
duce our method for creating the random potential and 
for simulating motion in this potential. Our results are 
presented in section 4. A discussion of the results is pre- 
sented in section 5, where a comparison is made with 
the renormalization group predictions. We conclude in 
section 6. 



2. THE QUENCHED DISORDER 

We consider the motion of one charged particle in a sea 
of quenched charges in two dimensions. The statistics of 
the motion of the particle is completely determined by 
the statistics of the quenched potential field that the par- 
ticle encounters. The quenched charges, which obey bulk 
charge neutrality, give rise to a charge-charge correlation 
function that vanishes as k 2 for small k in Fourier space. 
The potential, which is the convolution of the charge dis- 
tribution with the logarithmic Coulomb law, gives rise to 
a potential-potential correlation function that diverges as 
1/k 2 for small k |2^]. This long-ranged correlation func- 
tion for the potential experienced by the diffusing particle 
is exactly of the form that leads to anomalous diffusion. 

In two dimensions the interaction between the charges 
and the diffusing particle, and between the quenched 
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charges themselves, is logarithmic. Electrostatic charges 
in two dimensions, or line charges in three dimensions, 
interact with this law. Certain topological defects in two 
dimensions also interact with this same law. For exam- 
ple, dislocations in solids interact with a logarithmic law 
due to induced long-range elastic strain fields. Disclina- 
tions in hexatic membranes also interact with this effec- 
tive law, due to screening of the induced strain fields by 
free dislocations. Typical examples of these topological 
defects include line defects in three dimensional crystals, 
vortices in superfluids, and flux lines in superconductors 
11- 

The exponent for the mean square displacement de- 
pends indirectly on the density of defects via the prefac- 
tor of the potential-potential correlation function. The 
form of correlation function appropriate for small k is 
Xw(k) ~ l/k 2 , where 7 is the strength of the disorder, 
and X OT (k) = J d d r exp(zk • r)x„„(r) is the Fourier trans- 
form of the correlation function. We are free to chose dif- 
ferent behavior away from the origin in Fourier space. A 
natural choice for this correlation function is the inverse 
of the diffusive Green's function. On a square lattice with 
spacing Ar we, thus, use the form 



1 f d d k 



7(Ar) 



4 - 2cos(fe r Ar) - 2cos(fc y Ar) 



(1) 



Note that this form of Xw(k) has the appropriate lim- 
iting behavior as the lattice spacing become infinitesi- 
mally small and as the wavelength becomes large. The 
particles move through this potential field starting from 
random initial positions. Renormalization group calcu- 
lations have rather convincingly shown that the mean 
square displacement exhibits an anomalous behavior at 
long times |l]-||,^6|] . On a log- log scale, the mean square 
displacement as a function of time has a slope of 1 — 5, 
where 



5 = 



1 



8tt 



(2) 



and (3 = l/(k B T). 



3. SIMULATION METHOD 

We now consider the creation of the Gaussian random 
potential V(r) on a lattice. The potential takes on real 
values at each lattice site. The probability of observing 
any specific potential distribution is given by 



P[V] 



-0H[V] 



(3) 



where 



(3H[V] = \J d d vd d v'V(v)x^(v-v')V{v') 



2 J (^IWfcM 
= ^E^£l^ k )| 2 ^(k), (4) 

Z is a normalizing constant, and Afc = 2w/(NAr). Here 
the lattice is in d dimensions (d — 2 in our case), has 
iV" unit cells on a side, and has lattice spacing Ar. The 
Fourier transform is given by V(k) = J d d rV(r) exp(zk • 
r) = £ r (Ar)'V(r)exp(ik-r). Since V(-k) = V*(k), we 
have 



VL/2 



Afc 
2^T 



\V(k)\ 2 x - v \k) (5) 



with ki/2 meaning half of fc space and 

1/2, if k is on a special point 

-k) { , (6) 

1, otherwise 

where the special points are the origin, the corners of 
the lattice, and the intersections of fc^-axis and fc y -axis 
with the lattice boundaries. In this form, V^(k) and V"(k') 
are independent as long as k 7^ k'. The potential V(k) 
is Gaussian with the following variance for the real and 
imaginary components for k not special 



ReF(k) 



ImV(k) 



1 



(7) 



(8) 



where ft — (NAr) d is the volume of the system. If k is 
special, 



and 



JmV(k) = , 



ReV(k) )=nx w (k) 



(9) 



(10) 



Also, by definition, V(0) = J d d rV(r). Since it is not the 
actual magnitude of the potential but rather its gradient 
that is of interest, we define V — V — (V). In this form, 
V'(0) = and W = VV. We use the potential V 
in the simulation. To create the potential V'(r) that we 
need, we first create ^'(k) in half of k-space by generating 
independent Gaussian random numbers with variances 
given by Eqs. (|t|)-(^0|). We then generate the other half 
of k-space using the relation V'(-'k) = V'*(k). Finally, 
we generate V' in real space by performing an inverse 
fast Fourier transform fl27|]. 
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Alternatively, we can construct a real potential by first 
generating a complex V(r) and then extracting a real po- 
tential, V'(r). Specifically, we generate a complex Gaus- 
sian field with the probability distribution 



P[V] 



Z -1 exp 



J d d r<fW*(r)x~ 1 (r-r')V(r / ) 



(11) 



Note that the fields ReV(r) and lmV(r) are each Gaus- 
sian. Thus, we find 



We define 



(V*(r)V(r'))= X (r-r'). 



V'(r) = ReV(r) +ImV(r). 



(12) 



(13) 



V(r) is also Gaussian, since 



j d d rd d r'V*(r)xvv(r - r')V(r') 
= J d d rd d r'ReV(r) x -i (r - r')ReF(r') 
+ / d^rdVlmV^x^Cr-rOlmVfr') 
-i / d d rd d r'lmV{r) X yy(r - r')ReV(r') 
+i J d d rd d r'ReV(r) X ^(r - r')ImV(r'). (14) 

The complex part of the above equation vanishes, since 
Xvv( r - r ') = Xvv( r> - r ) for a medium obeying Eq. (Q). 
From our definition of V , we find 



(y'(r)F'(r'))= Xw (r-r'). 



(15) 



This definition of V leads to a real, Gaussian potential 
field with the correct statistics. It is, therefore, an equally 
valid quenched random potential. 

With the potential established, we shift focus to the 
initial conditions for the particle and the mechanics for 
diffusion. In the theoretical treatment of anomalous dif- 
fusion, the particles are assumed to be uniformly dis- 
tributed over the lattice. In our simulation, therefore, we 
choose the initial positions of the random walkers from 
a uniform distribution. In addition, we also examine the 
case of a Boltzmann distribution of initial positions, as 
would be appropriate, for example, for a NMR experi- 
ment on an equilibrated system. The particles are se- 
lected from a random initial position on a finite lattice, 
and they begin to diffuse at t = 0. We could let one parti- 
cle diffuse an infinitely long time and record its behavior 
exactly as prescribed by theory. Depending on whether 
the diffusion is self-averaging, we may also need to aver- 
age over different realizations of the potential. Unfortu- 
nately, lattice effects will appear at long times as a result 
of the periodic boundary conditions. Another method, 
which is more efficient and yields better statistics, is to 
sample many random walkers for a shorter period of time. 
This will work as long as the observation time is long 



enough to be in the scaling region, that is, as long as the 
system is large enough. We have investigated finite size 
effects and found that a square lattice of N = 2048 unit 
cells on a side is sufficiently large to allow the particles 
to be in the scaling regime for our range of parameter 
values. At short times, the particles will display normal 
diffusion behavior. At intermediate times, the diffusiv- 
ity will tend to zero. At long times, we will observe the 
anomalous sub-diffusion, where {r 2 (t)) is expected to be 
proportional to i 1_<5 as t — ► oo. 

The statistics of the random walk on the lattice are 
conveniently described by a master equation. This mas- 
ter equation defines the probabilities and rates of all pos- 
sible hops that the random walker can execute. We ex- 
actly solve this master equation by a Poisson process Q . 
In one dimension, the master equation looks like 



dP n 

dt 



U n -\P n -i — D n P n + D n+ iP n+ i — U n P n , (16) 



where P n is the probability for being at site x n at time 
t, U n is the rate at which transitions occur from x n to 
x n +i and D n is the rate at which transitions occur from 
x n to x n -\. In a two dimensional space, we denote the 
rates by T n ^ m , T n ^ m , T njTn , and T n<m , where T n<m and T n ,m 
are rates at which transitions occur from x n ^ m to Xn+i.m 
and x n -i tm respectively, and Tn)m and r^m are rates at 
which transitions occur from x nim to x n-m+ i and x n>m -\ 
respectively. The two dimensional master equation is 

dP n ,m _ (1) 



dt 



T W P . _ T (i) p 

' n-l,m r n-l,rn 'n,m r n,i 

+ r (2) P . i - P 

+ T ( 3 ) P . _ T 0) p 
' 'n,m— l r n,m— 1 'n,m r n,r 

+ T (4) P . i - T< 4 ) P 



(17) 



We demand that this master equation lead to a Boltz- 
mann distribution of the random walkers at long times. 
In other words, we want the stationary probability P^ of 
being at site x n ^ rn to be 



P 



(18) 



This distribution will arise if detailed balance is enforced. 
For example, the condition of detailed balance for tran- 
sitions in the x direction is 



-(2) 



n+l,m n+l,m — T n,m r n,m 



(1) P 

n on J n. 



(19) 



This criterion implies, with the use of the Boltzmann 
distribution, 



r n-\-l,m >n,m 

1 n.m x ' 

n+l,m 



-lPV(x n+ i, m )-V(x n , m )] 



(20) 
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Hence, one consistent expression for the rates Tn]m and 

(2) . 
7"n,m IS 



,l3[V(x n , m )-V(x n+1 , m )]/2 



(1) = D 
T n,m (Ar) 2 L 

(2) = D (3\V(x n . m )-V(x,^ 1 . m )]/2 



(Ar) 



(21) 
(22) 



Similar relations hold for the transition rates in the y 
direction: 



(3) = D r 0\V(x„. m )-V(x n , m+ i)]/2 
n, m ' * ^ " ° 



-(4) 



(Ar) 2 

£> 
(AO 



,/3[^n,m)-V(x„, m -i)]/2 



(23) 
(24) 



We use these transition rates to generate a stochas- 
tic Poisson process. Specifically, a particle begins at a 
site Xn ym . The particle waits at this site with an expo- 
nentially distributed amount of time characterized by its 
mean value 



(dt) 



1 



Jl) , (2) , (3) , (4) 
'n.TTi 'n.m 'n.m 'n.m 



(25) 



We generate the actual time increment via 

1 



T (D , (2) , (3) , (4) 
•n.m 'n.m \ 'n.m \ >n.m 



ln(x) , (26) 



where x is uniformly distributed random number with 
< x < 1. We then generate a second uniformly dis- 
tributed random number that we use to pick one of the 
four possible nearest-neighbor hops according to their 
probabilities: 







P(j^n,7n * ^n.m — 1 





r (1) 




T (1) - 
' 77, m 


hr (2) +r (3) 

1 / 77,777 \ ' 71,777 


+ T (4) 

T '71,777 




r (2) 

' 77,771 




" r (1) - 


hr (2) +r (3) 

r / 77,777 ' 71,777 


+ T (4) 

1 '77,777 




T (3) 

' 77 , 777 




" T (1) - 
' 77, m 


Vt {2) +r (3) 

1 / 77,777 \ 1 77,777 


+ T (4) 

1 '71,777 




_(4) 

' 77,777 




r (1) - 


hr (2) +r (3) 

1 / 77 , 777 ~ '71,777 


+ T (4) 

1 ' 77,777 



(27) 



The particle hops to the new site, and time is incremented 
by dt. 

On a finite lattice with periodic boundary conditions, 
particles reenter the lattice when they reach the bound- 
aries. Furthermore, the shortest path between two parti- 
cles may cross the boundary. This is significant, because 
we are fundamentally interested in the implications of 
anomalous diffusion for chemical reaction, and the ap- 
propriate measure for defining distance between two re- 
actants is the length of the shortest path. We, therefore, 
define the distance traveled by a diffusing particle as 



r 2 = min { [i - i Q + pN] 2 + [j - j + qN] 2 } (Ar) 2 . (28) 

p,q I > 

Here io and jo are the initial position coordinates of the 
particle, N is the dimension of the lattice, and the min- 
imum over integer p and q mathematically defines the 
shortest path. 

An important component of the simulation is the ran- 
dom number generator. A desirable generator ensures 
that the correct random statistics are used in lattice cre- 
ation, choice of initial particle positions, choice of hop- 
ing directions, and generation of time increments. In 
this study, we used two random number generators: one 
that is a sum of three linear congruential generators 0] 
and an exclusive-or, linear-feedback shift-register method 
fjOf . Unless specified otherwise, all the data that are 
shown were generated using the sum of three linear con- 
gruential generator method. 



4. RESULTS 

We are interested in the long time diffusive behavior of 
the particles diffusing on the disordered lattice. Since the 
diffusion is anomalous, the slope of ln(r 2 ) versus hit will 
not be unity. This slope will approach as the strength of 
disorder goes to oo and will approach 1 as the strength 
of disorder goes to 0. We, therefore, examine disorder 
strengths of varying magnitudes. We find that values of 
/3 2 7 in the range of 1 to 20 produce convincing and high 
quality results. For /3 2r y too small, the slope will be 1 — e, 
with e smaller than the noise in our simulation, making 
the anomalous diffusion difficult to observe. When /3 2 7 
is too large, strong fluctuations in the particle behavior 
appear at short times because of significant lattice effects. 
Within the chosen range of /3 2 7, we observe well behaved 
curves that exhibit distinct disorder strength dependent 
slopes. 

We perform several simulations with different initial 
starting positions for the random walker. We collect 
these data as histograms of (r 2 (i)) versus t, with a tem- 
poral bin width of to = 1. When plotted on a logarith- 
mic scale, these histograms will have much more data for 
large \n(t/to) than for small ln(t/to). In order to coun- 
teract this effect, we use an exponential sequence for se- 
lecting data from the histogram when performing the fit 
to determine the slope. Data for both short times, when 
the diffusion is not yet anomalous, and long times, when 
finite size effects are significant, was not used in the fit- 
ting procedure. This procedure leads to widely spaced, 
independent data points. A convenient byproduct of this 
procedure is that the standard error of the fit gives a 
reliable estimate of the error in the measured 8. These 
error bars are included in all of the figures. Note that 
there could also be a systematic error in each simula- 
tion related to the fact that only a single realization of 
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FIG. 1. Shown are the Monte Carlo results for the slope 
of the mean square displacement, m, as a function of strength 
of disorder, /3 2 7- We have used Eq. (^) to generate the lattice 
and we placed the random walker uniformly and randomly at 
the beginning of each random walk. Shown in dashed are the 
renormalization group predictions. 
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FIG. 2. The same quantities as in Figure 1, but using Eq. 
([l3]) to generate the random potential lattice. 
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the disorder is employed. This systematic error seems to 
be small, as simulations with different realizations of the 
disorder lead to values of 5 that differ approximately by 
the standard error of the fit. 

We pick enough different starting positions and follow 
each simulation long enough to ensure adequate statis- 
tics. We found that Njo = 10000 particles is sufficient 
to produce fairly smooth histograms for the mean square 
displacement. We also found that Ni en — 2000000 steps 
in each random walk is enough for the particles to reach 
the asymptotic scaling regime. We found that a lattice 
of size N = 2048 is sufficient to give results that span a 
broad range of times in the scaling regime. In all simula- 
tions, the lattice spacing Ar is set equal to unity, which 
we can enforce by a spatial rescaling. We also set the 
bare diffusivity equal to unity, which we can similarly 
enforce by a temporal rescaling. These rescalings will 
not affect the scaling behavior at long times. The value 
of S observed at long times, for example, is independent 
of these bare values. 

Figure 1 shows the slopes determined from simulations 
with strengths of disorder 1 < /3 2 7 < 20. We have used 
Eq. (|5|) to create the lattice in these simulations. Wc 
have chosen the initial position of the random walker uni- 
formly on the lattice, in direct correspondence to the case 
considered by the renormalization group studies. The 
data points represent slopes of the mean square travel 
distances of the random walkers as a function of time. 
We calculated these values by fitting \og(r 2 (t) /(Ar) 2 ) as 
a function of ln(i/io), as described above. For each value 
of the disorder strength, three independent runs were per- 
formed on different realizations of the disorder. We show 
in Figure 1 the average slope (circles) and the associated 
standard error of the fit (error bars) for each strength 
of disorder. The renormalization group predictions are 
shown as the dashed line. We see excellent agreement 
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FIG. 3. The same quantities as in Figure 1, except that 
we have placed the random walker randomly according to a 
Boltzmann distribution at the beginning of each random walk. 



between the simulation results and the theoretical pre- 
dictions, with nearly all the observed values within one 
standard deviation of the expected value. The varying 
standard deviations reflect randomness in the potential 
fields, initial positions, and hopping rates. The slopes 
do not include contributions from short time behavior 
or long distance behavior [(r 2 ) on the order of (TVAr) 2 ]. 
Both of these regimes exhibit significant lattice effects, 
effects not considered in the renormalization group stud- 
ies. 

Figure 2 shows similar data derived from lattice po- 
tentials constructed from complex fields, Eq. ([l3|). These 
results should be identical to those of Figure 1. We ob- 
serve that the average values are, again, consistent with 
the theoretical predictions. Interestingly, the standard 
errors are consistently smaller than those of Figure 1. 

Figure 3 shows the slopes that result when the particle 
initial positions are chosen from a Boltzmann distribu- 
tion. These conditions mimic those of a transient ex- 
periment, such as pulsed field gradient NMR, performed 
on an equilibrated system. We see agreement between 
the observed values and the renormalization group pre- 
dictions, indicating that this change in initial conditions 
is not "relevant" in the technical sense. We do see less 
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FIG. 4. The same quantities as in Figure 1, but using the 
the correlation function Xw(k) = 7exp( — k 2 /2)/k 2 in creat- 
ing the random potential. 

consistency in the standard deviations due to the non- 
uniform sampling of the rugged potential landscape. 

Figure 4 shows the slopes that result when one uses 
the correlation function 

-k 2 II 

X^=1^-- (29) 

This correlation function has the same small k behav- 
ior as Eq. (|l|), but distinct behavior for large k. The 
long-time behavior of these two correlation functions are 
expected to be the same, as long as /3 2 7 is not renormal- 
ized by changes in the large k behavior. The prefactor of 
the mean square displacement, however, is observed to 
be significantly larger for Eq. ( p9| ) than for Eq. (Q). 

5. DISCUSSION 

The results of the computer simulations agree well with 
the predictions of the renormalization group studies. For 
each value of /3 2 7, the simulations yielded a slope of 
ln(r 2 /(Ar) 2 } versus \n(t/to) in excellent agreement with 
the analytical predictions. An important observation is 
that even at long times, /3 2 7 does not become renormal- 
ized. This is a non-trivial observation, as details in the 
simulation that are technically "irrelevant" could renor- 
malize f3 2 j a finite amount. For all values of /3 2 7, lattice 
effects produce normal diffusive behavior at short times. 
This normal diffusion crosses over to anomalous diffusion 
fairly quickly, within a time corresponding to relatively 
few hops by the particle. At very long times, the mean 
square displacement reaches a maximum value due to 
the periodic boundary conditions. This maximum value 
is proportional to iV 2 , since the mean square displace- 
ment defined in Eq. (^8|) is always less than or equal to 
(NAr) 2 /2. Figure | shows the mean square displace- 
ment measured in a typical run. The plateau at long 
times is clearly visible. The normal diffusion at short 
times crosses over to anomalous diffusion so rapidly that 
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0.0 5.0 10.0 15.0 

ln(t/t ) 

FIG. 5. Shown is the mean square displacement versus 
time for f3 2, y = 1 and Xw(k) = 7exp( — k 2 /2)/k 2 . We have 
used Eq. (|^) to generate the lattice and we placed the random 
walker uniformly and randomly at the beginning of each ran- 
dom walk. Note the plateau in the mean square displacement 
at long times. 

it is not visible with the histogram bin width we used 
(to = 1). 

As discussed, both methods for generating the ran- 
dom lattice potentials, Eq. (||) and (|l3|), lead to Gaus- 
sian random fields with the correct correlation function. 
Both methods should produce the same results for the 
mean square displacement. We find that, indeed, both 
approaches give the same average result for the mean 
square exponent. The real lattice generation is more effi- 
cient in terms of memory utilization. One implication of 
this efficiency, however, is that the real lattice generation 
is constructed from fewer independent random numbers 
and is a more severe test of the pseudo-random number 
generator. The fewer degrees of freedom used when im- 
plementing Eq. (|^) when compared to Eq. ([l3]) is most 
likely what leads to the larger error bars in Figure 1 when 
compared to Figure 2. 

The long-range correlations in the random potential al- 
low, in principle, for the possibility that the mean-square 
displacement depends on the distribution of initial condi- 
tions. All of the theoretical predictions, for example, are 
based upon the assumption of a uniform distribution of 
initial conditions. In experiments upon equilibrated sys- 
tems, however, the initial conditions are distributed in a 
Boltzmann manner. We see from Figure 3 that Boltz- 
mann initial conditions lead to the same mean square 
exponent at long times. 

Two forms of the potential-potential correlation func- 
tion, Eq. (|l]) and (p9|), are used to explore the depen- 
dence on irrelevant, large k details. At large distances 
and long times, the observed scaling behavior that re- 
sults from the two correlation functions would differ only 
if technically irrelevant details of the lattice renormalizc 
(3 2 j. We found that, indeed, the prefactor of the relation 
{r 2 (t)) ~ at 1 ^ 8 does depend on these irrelevant details. 
As we see from Figure 4, however, the exponent is inde- 
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FIG. 6. Shown is the mean square displacement versus 
time for /3 2 7 = 10 and Xvv(k) = 7 exp( — k 2 /2)/k 2 . We have 
used Eq. (H) to generate the lattice and we placed the ran- 
dom walker uniformly and randomly at the beginning of each 
random walk. The top solid curve comes from using the ex- 
clusive-or, linear-feedback shift-register pseudo-random num- 
ber generator and the bottom solid curve comes from 
the sum of three linear congruential pseudo-random number 
generators pjl. The bottom curve has the expected slope of 
0.71 = 1 - 1/(1 + 871-/IO) at long times. The top curve has 
a slope of 1.09 at the longest times shown, which indicates 
superdiffusion. Dashed lines with slopes of 0.71 and 1.09 are 
shown for convenience. 



pendent of these details. 

The choice of pseudo-random number generator to use 
is an important consideration in all Monte Carlo sim- 
ulations. A pseudo-random number generator is used 
in three components of the present simulation: in lat- 
tice creation, in selecting the particle initial positions, 
and in creating the transition rates and hop directions. 
The two random number generators employed in our sim- 
ulation are a sum of three linear congruential genera- 
tors pj| and an exclusive-or, linear- feedback-shift reg- 
ister method with table length 55 |3(]]. Both of these 
generators are thought to be fairly reliable. All of the 
results in Figures 1-5 were produced by the sum of three 
linear congruential generators method. Figure 6 com- 
pares the mean square displacements produced by these 
two generators. The three linear congruential genera- 
tor method results in slopes for all values of f3 2 j that 
are consistent with theory. Interestingly, the linear feed- 
back shift register method always leads to slopes greater 
than expected. Using this generator, one would conclude 
that j3 7 is renormalized a finite amount, by some un- 
known factor. In fact, the factor is the inadequacy of the 
pseudo-random number generator! At long times, the ex- 
ponent of the mean square displacement exceeds unity. 
This super-diffusive behavior is in conflict with a rigor- 
ous bound known for diffusion in Gaussian random me- 
dia: Ym H ^{r 2 {t)/t) < 4Dexp[-/? 2 XTO (0)] @. In fact, 
at long times, this linear feedback shift register method 
appears to produce ballistic behavior, (r 2 (t)) ~ at 2 (data 
not shown). This incorrect super-diffusive behavior ap- 



pears only for potentials that lead to anomalous diffu- 
sion. The long-range correlations in the potential, which 
lead to the anomalous diffusion, apparently couple to the 
residual correlations in this pseudo-random number gen- 
erator. Similar, poor results from this type of generator 
have been observed in another system with long-range 
correlations — the Ising model at its critical point f$^| . In 
this feedback generator with a short table length 

led to predicted critical exponents differing from the true 
values by many times the estimated standard deviation. 



6. CONCLUSION 

We find a satisfying match between theory and nu- 
merical results in our simulation of anomalous diffusion. 
We found that anomalous behavior occurs in the long 
time regime, with a transition from normal diffusion at 
short times. We find that the prefactor for the scaling of 
the mean square displacement is renormalized by short- 
distance correlations in the potential, although the expo- 
nent is not. Reasonable distributions for the initial con- 
ditions lead to the same exponent for the mean square 
displacement. Interestingly, the correct anomalous diffu- 
sion behavior is observed only with a high quality pseudo- 
random number generator. 
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